// This file is part of Spectrum Editor Called SECS.
//
// Spectrum Editor Called SECS is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// Spectrum Editor Called SECS is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with Spectrum Editor Called SECS.  If not, see <http://www.gnu.org/licenses/>.

typedef std::complex<double> Complex;

struct Signal {
  int sample_rate;
  std::vector<Complex> data;

  Signal() { }

  Signal(int _sample_rate) {
    sample_rate = _sample_rate;
  }

  Signal(int _sample_rate, std::vector<Complex> _data) {
    sample_rate = _sample_rate;
    data = _data;
  }

  Signal(Wave wave) {
    sample_rate = wave.sample_rate;
    for (int i = 0; i < wave.samples.size(); i++)
      data.push_back((float)wave.samples[i][0] / wave.maximum_amplitude);
  }

  void _fft(std::vector<Complex>& _data) {
    int N = _data.size();
    if (N < 2)
      return;
    std::vector<Complex> even, odd;
    for (int i = 0; i < N; i ++) {
      if (i % 2)
       odd.push_back(_data[i]);
      else
        even.push_back(_data[i]);
    }
    _fft(even);
    _fft(odd);
    for (int i = 0; i < N / 2; i++) {
      Complex t = std::polar(1.0, -2.0f * M_PI * i / N) * odd[i];
      _data[i] = even[i] + t;
      _data[i + N / 2] = even[i] - t;
    }
  }

  Signal fft(int window_size, int fft_size) {
    std::vector<Complex> _data (data.begin(), data.begin() + window_size);
    for (int i = 0; i < fft_size - window_size; i++)
      _data.push_back(0);
    _fft(_data);
    return Signal(sample_rate, std::vector<Complex> (_data.begin(), _data.begin() + fft_size / 2));
  }

  // TODO ifft - how to handle zero padding?

  Signal to_polar() {
    Signal new_signal(sample_rate);
    for (int i = 0; i < data.size(); i++)
      new_signal.data.push_back(Complex (abs(data.at(i)), arg(data.at(i)))); 
    return new_signal;
  }

  // TODO to_cartesian

  Signal window(std::vector<double> window) {
    Signal new_signal(sample_rate);
    for (int i = 0; i < window.size(); i++)
      new_signal.data.push_back(data[i] * window[i]);
    return new_signal;
  }

  Signal get_segment(int begin, int end) {
    std::vector<Complex> _data(data.begin() + begin, data.begin() + end);
    return Signal(sample_rate, _data);
  }

  double root_mean_square() {
    double sum_of_squares = 0.0;
    for (int i = 0; i < data.size(); i++)
      sum_of_squares += real(data.at(i)) * real(data.at(i));
    return sqrt(sum_of_squares / data.size());
  }

  double threshold_curve(double frequency) {
    double b = 0.0075;
    double a_l = -24.0;
    double a_r = 32.0;
    return a_l - a_r + a_r * pow(b, frequency / 20.0);
  }

  int find_maximum() {
    int maximum = 1;
    for (int i = 2; i < data.size(); i++) {
      if (real(data.at(i)) > real(data.at(maximum))) {
        maximum = i;
      }      
    }
    // TODO quadratic interpolation
   return maximum;
  }

  void print() {
    for (int i = data.size() - 1; i > -1; i--)
      std::cout << data.at(i) << std::endl;
  }

  void histogram(int width) {
    for (int i = data.size() - 1; i > -1; i--) {
      int number_of_characters = (int)(real(data[i]) * width);
      for (int j = 0; j < number_of_characters; j++)
        std::cout << "*";
      std::cout << std::endl;
    } 
  }

  Signal normalize() {
    Signal new_signal(sample_rate);
    double maximum = abs(data.at(0));
    for (int i = 1; i < data.size(); i++) {
      if (abs(data.at(i)) > maximum)
        maximum = abs(data.at(i));
    }
    for (int i = 0; i < data.size(); i++)
      new_signal.data.push_back(data.at(i) / maximum);
    return new_signal;
  }

  Signal amplitude_to_decibel(double reference) {
    Signal new_signal;
    for (int i = 0; i < data.size(); i++)
      new_signal.data.push_back(Complex(10.0 * log10(real(data[i]) / reference), imag(data[i])));
    return new_signal;
  }
  
  void subtract_signal(Signal other_signal) {
    if (data.size() == other_signal.data.size()) {
      for (int i = 0; i < data.size(); i++)
        data[i] -= other_signal.data[i];
    }
  }

  Complex operator[](int i) {
    return data.at(i);
  }
};

